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Abstract. Numerical methods are proposed for an advanced Poisson-Nernst- 
Planck-Fermi (PNPF) model for studying ion transport through biological 
ion channels. PNPF contains many more correlations than most models and 
simulations of channels, because it includes water and calculates dielectric 
properties consistently as outputs. This model accounts for the steric effect 
of ions and water molecules with different sizes and interstitial voids, the 
correlation effect of crowded ions with different valences, and the screening 
effect of polarized water molecules in an inhomogeneous aqueous electrolyte. 
The steric energy is shown to be comparable to the electrical energy under 
physiological conditions, demonstrating the crucial role of the excluded vol¬ 
ume of particles and the voids in the natural function of channel proteins. 
Water is shown to play a critical role in both correlation and steric effects 
in the model. We extend the classical Scharfetter-Gummel (SG) method for 
semiconductor devices to include the steric potential for ion channels, which 
is a fundamental physical property not present in semiconductors. Together 
with a simplified matched interface and boundary (SMIB) method for treat¬ 
ing molecular surfaces and singular charges of channel proteins, the extended 
SG method is shown to exhibit important features in flow simulations such as 
optimal convergence, efficient nonlinear iterations, and physical conservation. 
The generalized SG stability condition shows why the standard discretization 
(without SG exponential fitting) of NP equations may fail and that divalent 
Ga^’*' may cause more unstable discrete Ga^"*" fluxes than that of monovalent 
Na"*". Two different methods — called the SMIB and multiscale methods — are 
proposed for two different types of channels, namely, the gramicidin A channel 
and an L-type calcium channel, depending on whether water is allowed to pass 
through the channel. Numerical methods are first validated with constructed 
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models whose exact solutions are known. The experimental data of both chan¬ 
nels are then used to verify and explain novel features of PNPF as compared 
with previous PNP models. The PNPF currents are in accord with the exper¬ 
imental I-V (V for applied voltages) data of the gramicidin A channel and I-C 
(C for bath concentrations) data of the calcium channel with 10“®-fold bath 
concentrations that pose severe challenges in theoretical simulations. 


PACS number(s): 87.10.Ed 02.70.-c 47.61.Cb 05.20.Jj 


I. INTRODUCTION 


The literature on numerical methods for drift-diffusion (DD) or Poisson-Nernst- 
Planck (PNP) models of semiconductor devices and ion channels is large, in¬ 
cluding [T1l21l5|TfS1l6H7|5|91ITU|ll|ll2|ll3|ll4|ll5|ll6|ll7|ll8lll9ll20|l21 
therein. In biological simulations, continuum models have been challenged 


and references 
een challenged 

as inaccurate compared to Monte Carlo (MC), Brownian dynamics (BD), or 
molecular dynamics (MD) due to the gross approximation of atomic properties 
of channel proteins and electrolyte solutions [231I24II25II26II27II28II29II30II31II32II33] . 
Continuum models on the other hand have substantial advantages in efficiency 
that are of great importance in studying a range of conditions and concen¬ 
trations especially for large nonequilibrium or inhomogeneous systems, as are 
present in experiments and in life itself 


Based on the configurational entropy model |12] for aqueous electrolytes with 
arbitrary K species of nonuniform size, hard spherical ions, we extended the 
Poisson-Fermi model in [12] to a new model — called the Poisson-Nernst- 
Planck-Fermi (PNPF) model — for nonequilibrium systems by including specif¬ 
ically the excluded volume effects of the next species (K+l) of water molecules 
and the interstitial voids {K + 2) between all particles [13]. The PNPF model 
differs from most channel models in several respects: (i) it computes dielectric 
properties as an output that in fact vary with position and with experimental 
condition, (ii) a fourth order Cahn-Hilliard type partial differential equation 
emerges to replace the second order Poisson equation of PNP, which has a 
richness of behavior beyond the usual second order PNP description, and (iii) 
using the methods of this paper, this more powerfully correlated model is 
in fact much easier to compute in three dimensions than other steric PNP 
models. Previous work [43] gives more details. 


The PNPF model also provides a quantitative mean-field description of the 
charge/space competition mechanism of particles within the highly charged 
and crowded channel pore. The steric energy lumps the effects of excluded 
volumes of all ions, water, and voids. It yields an energy landscape of ions 
that varies significantly with bath concentrations in a 10®-fold range of exper¬ 
imental conditions for L-type calcium channels. 
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Computational Challenges. The 10®-fold range of bath conditions and the 
highly energetic behavior of permeating ions throngh the extremely crowded 
narrow channels pose severe challenges in implementations. The strength of 
local electric fields in a calcinm channel can be higher than that in a semi- 
condnctor device (comparing, for example, 0.27 V/nm estimated from Fig. 13 
in this paper for ion channel and 0.06 V/nm from Fig. 7 in |T8] for semicon- 
dnctor device). This means that the convergence and stability problems in ion 
channel simnlation can be more severe than those in semicondnctor devices. 
These problems are not made easier by the presence of conntervailing steric 
potentials of the same order of magnitnde. 

Moreover, geometric complexity and singnlarities of molecnlar snrfaces sepa¬ 
rating electrolyte solntions from protein atoms in biological systems need to 
be carefnlly treated in order to obtain tolerable accnracy in 3D PNP simnla- 
tions [22]. Seemingly small nnmerical over approximations can lead to errors 
that make resnlts not nsefnl. A second-order method called the matched in¬ 
terface and bonndary (MIB) method was developed by Wei et al. [22E1] for 
Poisson-Boltzmann and PNP models and is simplihed (SMIB) in [35] for the 
PF model to deal with the geometric singnlarities by the standard hnite dif¬ 
ference approximation. 

The Scharfetter-Gnmmel (SG) [2] method is an optimal and nniformly con¬ 
vergent method (with respect to the mesh size) to discretize drift-diffnsion (or 
Nernst-Planck) eqnations for flnx calcnlations becanse it integrates the corre¬ 
sponding ID initial valne problem exactly at every grid point [5]. We extend 
the classical SG method to the NPF eqnation by showing how the Fermi distri- 
bntion of hard spheres of water and ions is imposed. If the classical Boltzmann 
distribntion is nsed, the density of point charges wonld grossly overestimate 
ionic concentrations (that are in fact limited becanse of the hnite size of ions) 
and conseqnently lead to inaccnrate electrostatic potential and ion mobility 
by the classical PNP [2^291135] . We also show that the classical Goldman- 
Hodgkin-Katz hux approximation [36] in ion channels is in fact exactly the 
Scharfetter-Gnmmel hnx approximation on grid points in semicondnctor de¬ 
vices. Similar results appear in the seminal work of Mott m that was well 
known to Hodgkin, Gole, and Goldman. The pioneers in two different helds 
had the same idea that made a profound impact on their respective helds and 
others. 

The SG stability condition — a critical condition of the hux equation in im¬ 
plementation — is also extended to include the steric potential that is not 
present in classical PNP models. This stability condition explains why the 
standard hnite diherence or hnite element discretization fails when the elec¬ 
tric and/or steric potentials vary sharply in a layer region and the mesh of grid 
points is not sufficiently resolved. It plays a key role in preserving physically 
positive concentrations and divergence-free currents (current conservation) in 
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approximation [TO]. We take a closer look at the numerics concerning the ex¬ 
tended SG condition and discover that this condition is harder to satisfy for 
the standard methods for the divalent Ca^’*' than the monovalent Na"*" flux 
since the SG condition depends on the valence of ions. This is physically rea¬ 
sonable because Ga^^ ions are more energetic in binding and permeation in 
voltage-gated calcium channels that conduct Ga^"*" ions with high-hdelity and 
high-throughput [T8] . 

The combined method — the SMIB-SG method — is shown not only to achieve 
second order of accuracy for the PNPF model (with constructed exact solu¬ 
tions) but also to outperform the primitive SMIB method (without SG) for 
the gramicidin A channel due to the exactness property of the SG exponential 
htting between grid points. We also show that the primitive SMIB method 
fails to converge for the calcium channel due to its highly charged (-4e, e is 
the proton charge) and very narrow (about 1 A in radius) binding site as 
compared with that of the gramicidin A channel (-2e and about 2 A in ra¬ 
dius). In our simulations, water (1.4 A in radius) is found to flow through the 
gramicidin A channel but not to flow through the calcium channel in some 
conditions. We use a second method — called the multiscale method — that 
treats water and ions explicitly in the binding site of the calcium channel so 
that water may not move through the channel. It is multiscale since both 
Poisson’s theory of continuous charges and Goulomb’s law of discrete charges 
are used in the solvent domain. This demonstrates the novelty of the PNPF 
model as compared with previous PNP models in dealing with ion-protein, 
ion-ion, and ion-water interactions and the steric effect of ions and water in 
the narrow pore. PNPF captures many more of the correlations not present 
in PNP itself. It captures steric interactions of ions and water and packs them 
well (i.e. consistently) because it includes free space. Dielectric properties vary 
with position and concentration and are fully consistent with the rest of the 
model because they are outputs of the calcualtions, not inputs, as assumed in 
most channel models. 

The nonlinear algebraic systems of discrete PNP equations are very difficult to 
solve due to strong nonlinearity of the coupled system in both semiconductor 
devices and ion channels, especially with sharp potentials at practical applied 
voltages [5II?II1UII14II15II16II18II19II2UII21II22| . The PNPF model consists oi K + 2 
PDFs (1 fourth-order Poisson-Fermi and A' -|- 1 second-order Nernst-Planck). 
The fourth-order PF equation was proposed to account for the correlation 
effect of ions in water [39] and transformed to two second-order PDFs for 
computational efficiency and for calculating variable permittivity within the 
channel pore |42M5] . The last NP equation describes the dynamics of water 
molecules that play a critical role not only in the steric arrangement of all 
particles but also in its screening and polarization effects on ions in the sys¬ 
tem [MS] . The full PNPF model incorporates these atomic properties and 
thus can provide more accurate simulations but obviously at the expense of 
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more difficulties in implementation than that of previous PNP models. It is 
impractical to solve the {K + 2)M nonlinear algebraic equations resulting 
from a discretization of PNPF using Newton’s iteration on the coupled sys¬ 
tem, where the matrix size M corresponding to each PDE can easily grow 
to millions in 3D implementations. With a linearized Poisson equation, Gum- 
mel’s iteration is an efficient method because it solves each PDE successively 
[9]. It has been shown in |22] that an SOR-like method (without linearization) 
converges faster than Gummel’s method at higher bath concentrations for ion 
channel simulations provided that the relaxation parameter is appropriately 
chosen. We present a new SOR-like method for the PNPF model that differs 
from the previous models in the fourth-order PDE, the water NP equation, 
and the steric potential. It is shown that the method improves the convergence 
rate using the same gramicidin A channel protein as considered in |22] . 

The rest of the paper is organized as follows. In Section 2, we briefly describe 
the PNPF theory. Sections 3 and 4 present the numerical methods proposed 
in this paper. Section 4 also include two algorithms with respect to the SMIB 
and multiscale methods to illustrate implementation procedures for studying 
two different types of ion channels. In Section 5, the SMIB and SG methods 
are hrst validated by using a real protein structure of the GA channel with 
a set of exact solutions constructed for the PNP model. Both methods are 
shown to achieve optimal results as analyzed in this paper. The extended SG 
condition is then carefully scrutinized in discretization and used to explain 
why the standard discretization method is not feasible for the calcium channel 
model considered here, especially to approximate the high energetic Ga^^ flux. 
PNPF results are shown to agree with experimental I-V and I-G data of GA 
and calcium channels using the two algorithms, respectively. Some concluding 
remarks are made in Section 6. 

II. POISSON-NERNST-PLANCK-FERMI MODEL 


For an electrolyte in a solvent domain with arbitrary K species of ions and 
the next species iF -|-1 of water, the conhgurational entropy model proposed in 
[42] is extended in [43] to treat all particles as hard spheres with nonuniform 
sizes and to include explicitly as its last species K + 2 the voids between all 
particles. Based on the extended entropy model, the following Gibbs-Fermi 
free energy functional of the system is proposed in [43] 


^Fermi 


In (^V(r))^ - Y |V0(r)|^ + p(r)0(r) + 


g = ksT 


( K+l 

E 


V.=i 


C'i(r) In {vjCjij)) - Cjij) 

-Cj{r) In {vK+ 2 CK+ 2 {r)) - 


5 


( 1 ) 


where = e^co, e* is the dielectric constant of bulk water, cq is the vacuum 
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permittivity, Ic is a correlation length [191I50] . 0(r) is the electrostatic potential 
function of spatial variable r G p(r) = is the charge density, 

Cj{r) is the concentration of type j particles carrying the charge qj = ZjC with 
valence Zj and having the volume Vj = dvra^/S with radius a^, ks is the Boltz¬ 
mann constant, T is the absolute temperature, and /if = kBT\n{viCf/T^) 
is a constant chemical potential. Water is treated as polarizable spheres with 
zero net charge, so zk+i = Qk+i = 0. 


The total volume V of the system consists of the volumes of all particles and 
the total void volume vk+ 2 , i-e., V = VjNj +vk+ 2 , where Nj is the total 

number of type j particles. Under the bulk condition, dividing this equation 
by V yields the bath void volume fraction 


pB ^ '^K+2 

V 


K+1 


V 


K+l 

= 1 - E 

i=i 


( 2 ) 


where Cf is the bath concentration. The void fraction function 


r(r) 


1 - 


K+l 

E 


J=1 


^k+2C'k+2(j^), 


(3) 


varies with concentrations Cj(r) of all particles and thus with the distribution 
Cx+ 2 (r) of interstitial voids. 


Minimizing the Gibbs-Fermi functional (1) with respect to 0 and Gj yields the 
Poisson-Fermi equation [d2]|45ll49ll50] 

e, - l) V^(^(r) = ^ qiCi(r) = p(r) (4) 

i=l 

and the Fermi distribution 

C.(r) = Cfexp(-ft.#,(r) + S“(r)). S"(r) = In (5) 

respectively, where {3i = qi/(ksT) and is called the steric potential. The 

fourth-order PF equation reduces to the classical Poisson-Boltzmann (PB) 
equation —= p and the Fermi distribution reduces to the Boltzmann 
distribution Ci = Cf exp when Ic = = 0. The distribution (5) 

is of Fermi type since all concentration functions are bounded above, Gi(r) < 
1/vi [l3], i.e., Gi(r) cannot exceed the maximum value 1/vi for any arbitrary 
(or even inhnite) potential 0(r) at any location r in the domain 

If /c 7 ^ 0, the dielectric operator e = es(l —/^V^) approximates the permittivity 
of the bulk solvent and the linear response of correlated ions [50] . The dielectric 
function e(r) = €«/(! -|- rj/p) is a further approximation of e. It is found by 
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transforming (4) into two second-order PDEs [l5] 


( 6 ) 

(7) 


PFl : e, (/^V" - l) ^(r) = p(r) 

PF2 : VV(r) = ^(r) 

by introducing a density like variable 4/ that yields a polarization charge den¬ 
sity r] = —es4/ — p of water using Maxwell’s hrst equation |12]. Numerical 
approximation of the fourth order equation (4) was simplihed to the standard 
7-point hnite difference approximation of the second order equations (6) and 
(7) in [IS]- Boundary conditions of the new variable 4/ on the solvent bound¬ 
ary dQs were derived from the global charge neutrality condition [IS]- These 
functions make dielectric properties outputs in our model and calculations, 
unlike in most other treatments of channels. 

Including the electrostatic effect of a total of Q hxed atomic charges qj located 
at Tj in the biomolecular domain Qm that contains both channel protein and 
membrane lipids, the PF equation (4) is written as 


- l) VV(r) = II gi<5(r - r^) + H = p(r), Vr G 12, (8) 




Z=1 


where = 12^ U and 5(r — r^) is the delta function. Note that e = em^Q, 
Ic = 0, p(r) = Y .%1 - ^j) in and e = e^eo, h 7^ 0, p(r) = Y.f=i qiCi{r) 

in f2s, where is the dielectric constant of biomolecules. As mentioned above, 
numerical implementation of Fq. (8) (or Fqs. (6) and (7)) is complicated by the 
complex molecular surface dflm in real protein structures on which suitable in¬ 
terface conditions for the unknown functions T(r) and 0(r) should be properly 
imposed [13] . The approximation of interface conditions is not straightforward 
[ 2 ^ 441145 ] and can be made much worse by geometric singularities of dflm if 
the singularities are not properly treated. It was shown in [ 51 ] that the stan¬ 
dard, second-order hnite difference method is degraded to only by 

this kind of singularities, where h is the mesh size of grid points. 


For nonequilibrium systems, the classical Poisson-Nernst-Planck model [3^1133113^ 
can then be generalized to the Poisson-Nernst-Planck-Fermi model by coupling 
the hux density equation (in steady state) 


- V • Ji(r) = 0, r G f2s 


(9) 


of each particle species i = 1, ■ ■ ■ , K + 1 (including water) to the PF equation 
(8), where the hux density is dehned as 


J,(r) = -A VA(r) + Aa(r)V0(r) - A(r)VA’'^(r) 


( 10 ) 
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and Di is the diffusion coefficient. The flux equation (9) is called the Nernst- 
Planck-Fermi equation because the Fermi steric potential S'*’''^(r) is introduced 
to the classical NP equation. The NPF equation (9) reduces to the Fermi 
distribution (5) at equilibrium ^5] . 

The gradient of the steric potential in (10) represents an entropic force of 

vacancies exerted on particles. The negative sign in —means that the 
steric force is in the opposite direction to the ‘diffusionah force VCi, i.e., 

the larger = In (meaning more space available to the particle as im¬ 
plied by the numerator) at r in comparison with that of neighboring locations, 
the more the entropic force pushes the particle to the location r. The entropic 
force is simply opposite to the diffusional force VC* that pushes the particle 
away from r if the concentration is larger at r than that of neighboring loca¬ 
tions. Moreover, the Nernst-Einstein relationship |16] implies that the steric 
flux is greater if the particle is more mobile. Therefore, the gradi¬ 

ents of electric and steric potentials (V0 and VS'*’’'^) describe the charge/space 
competition mechanism of particles in a crowded region within a mean-field 
framework |13] . For more physical and mathematical details about the PNPF 
theory, we refer to |l3]. 

III. A GENERALIZED SCHARFETTER-GUMMEL METHOD 


We use the standard 7-point finite difference (FD) scheme in 3D [15] to dis¬ 
cretize the PNPF model. For ease of notation, we omit the subscript i in (9) 
when no confusion should arise. For conciseness, the FD discretization is sim¬ 
plified to ID in the following discussions as the corresponding 3D case follows 
obviously in a similar way. Furthermore, we only provide the FD formula for 
the flux equation (9) as the FD formulas with the SMIB method across the 
molecular surface dflm for Eqs. (6) and (7) have been given in [15], i.e., we 
consider 


dJ{x) 

dx 


A. 

dx 


-D{x) 


/dC{x) 
y dx 


+ /3C{x) 


d4>{x) 

dx 


C{x) 


dx J 


0 . ( 11 ) 


The primitive FD approximation of (11) is 

O'i-lCi-i + CliCi + Oj+iCj+i 
Ax'^ 

where Ax = Xj+i — Xi = h is the mesh size of a uniform grid on the x-axis in 
the domain, Ci ~ C{xi) is the unknown approximation of the concentration 










function C{x) at any grid point Xj, and the coefficients are given as 


a*-i = A-i [l - /5A0i_i/2 + A5 Ai/2 
= fli-1 + fli+1 ~ 2(Zi)j_l + A+i) 

a.+i = A+i [1 + /9A0,/2 - AA''72] , 


(13) 


where A^-i = <Pi — ~ ^i+\ — + a:i)/2 etc. The diffusion 

coefficient function D[x) is equal to a constant in the bath and to a reduced 
constant OD^ in the channel pore with 0 < 0 < 1. The function D{x) along the 
channel axis is constructed by using the interpolation method presented in [22] 
for connecting the bath value and the pore value OD^ such that D{x) is a 
continuously differentiable function. The factor 9 is the only tuning parameter 
in the PNPF model to £t experimental data |15ll22ll28ll35ll36ll40ll55ll56j . We shall 
investigate the magnitude of 9 for GA channel and compare it with those 
obtained by MD and BD simulations. The comparison is used to verify the 
correlation and steric effects considered in PNPF. 


At any two adjacent grid points x* and Xj+i, the FD approximation of the zero 
flux (J(x) = 0) is 


A+i - A ^ A+i + A , aA" ^ 

Ax 2 \ Ax Ax ) ’ 


(14) 


which implies that we may obtain the inequality 


A+i — Ci> A+i + Ci (15) 

and thereby a negative {unphysical) concentration A < 0 at x, if 

1 (-/ 3 A.#.. + ASp) > 1 . ( 16 ) 


Without the steric term AS'*’"'^, this inequality is the well-known Scharfetter- 
Gummel stability condition in semiconductor device simulations [2117] 

2 1^ 

- A(j)i = -(A+i -4>i) <- = —— (for /? > 0) (17) 

q 

required to ensure that the FD equation (12) does not produce unphysical 
approximations. Note that q = 2e for Ga^"*" yields an upper bound ksT/e in 
(17), which is a half of that for Na"*" and means that if the potential difference 
—A0j between two adjacent points is greater than ksT/e ~ 25.7 mV at room 
temperature, the resulting approximation of the Ga^"*" flux Jca2+ in (9) niay be 
completely unphysical, although the same discretization of the Na"*" flux JNa+ 
may still be feasible. In other words, the FD formula (12) is more unstable 
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for than for Na"*". In fact, if the SG condition is violated, Newton’s 

iteration for solving the coupled PNP system of nonlinear algebraic equations 
is generally divergent. Of course, we could reduce the mesh size /S.x so that the 
difference —A0j is small enough to satisfy (17) at all grid points i. This would 
however incur larger algebraic systems (and thus larger conditioning numbers 
of the system) for which the computational cost would be more expensive. 
Even using adaptive meshes that efficiently resolve internal or boundary layer 
regions where —A^j varies sharply, the primitive approximation (12) would 
still diverge or show extremely slow convergence if the layer thickness is very 


The convergence and stability issues are further complicated by the steric 
potential 5'*’’'’ in ion channel simulations if it is added to the FD flux equation 
(12) as given in (13). From (16), we obtain a new SG condition for ion channels 


- /3A0i + A^f ^ < 2 


( 18 ) 


that will be a focal point in our numerical investigations in Section 5. 


Stabilization. To stabilize (12), we extend the classical Scharfetter-Gummel 
approximation [2] of the flux J(x) to include the steric potential such that 



D 

Ax 


lB(-U)a+i - B(U)Q] 


( 19 ) 


where U = (5A(j)i — ASf‘^ and B{t) = is the Bernoulli function [7]. Eq. (19) 
is an exponential htting scheme for the concentration function C{x) between 
the mesh points Xi and Xj+i and is derived from the assumption that the flux 
J, the local electric held and the local steric held ^re all constant in 
this subinterval, i.e., 

J ^ dC{x) _ X e (xi, Xi+i), (20) 


where k = Solving this ordinary diherential equation (ODE) with a 

boundary condition G, or Gj+i yields the well-known Goldman-Hodgkin-Katz 
hux equation in ion channels |16], which is exactly the same as that in (19). 
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The generalized Scharfetter-Gummel method for (11) is thus 


dJ{Xi) ^ + QiCi + Oi+lCi+l _g /2j^\ 

dx Ax 

J,-i = ^ 

J.+i = ^ |B(-t.)Q+i - B((i)Q] 

(i = - AS*", B(() = Gt 

— 1 

ai-i = Qi = + B{ti), Oj+i = -B{-ti). 

The SG method is optimal in the sense that it integrates the ODE (20) ex¬ 
actly at every grid point with a suitable boundary condition [5]. Therefore, 
the SG method can resolve sharp layers very accurately [5] and hence needs 
few grid points to obtain tolerable approximations when compared with the 
primitive FD method. Moreover, the exact solution of (20) for the concentra¬ 
tion function C{x) yields an exact flux J(x). Gonsequently, the SG method is 
current-preserving, which is particularly important in nonequilibrium systems, 
where the current is possibly the most relevant physical property of interest 
It is difficult to overstate the importance of the current preserving fea¬ 
ture and it must be emphasized for workers coming from fluid mechanics that 
preserving current has a significance quite beyond the preserving of flux in 
uncharged systems. The electric field is so strong that the tiniest error in pre¬ 
serving current, i.e., the tiniest deviation from Maxwell’s equations, produces 
huge effects. The third paragraph of Feynman’s lectures on electrodynamics 
makes this point unforgettable |S7|. Thus, the consequences of a seemingly 
small error in preserving the flow of charge are dramatically larger than the 
consequences of the same error in preserving the flux of mass. 


IV. SMIB-SG AND MULTISCALE METHODS 


To test the PNPF theory and verify the numerical methods developed in 
this paper, we consider the GA channel with a real protein structure and a 
simplihed molecular model of L-type calcium channels. The main difference 
between these two channels is that the GA channel has a more rigid and less 
negatively charged pore with about 2 A in radius whereas the Ga^+ channel 
has a flexible and higher negatively charged binding site with radius varying 
from 1 A to 2.5 A. The GA channel is also much longer (22 A, see below) than 
the selectivity filter of the L type calcium channel (10 A [58] )• Gonsequently, 
the GA channel is only cation selective whereas the Ga^+ channel is exquisitely 
Ga^+ selective. The steric potential is a key component of PNPF to properly 
describe this important difference in selectivity along with the size effect of 
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(a) 


(b) 


Fig. 1. (Color online.) (a) Top view of the gramicidin A channel, (b) A cross section 
of 3D simulation domain for the GA channel. The channel is placed in a cubic box 
with the length of each side being 40 A and the thickness of the membrane being 
24 A. 

water (1.4 A in radius). We use two different treatments of water that yield 
two different steric potentials and size effects. 


A. The SMIB-SG Method for Gramicidin A Channel 


Fig. 1(a) is a top view of the GA channel downloaded from the Protein Data 
Bank [59]. A 2D cross section of the 3D simulation domain of the channel 
embedded in a membrane is sketched in Fig. 1(b), where the biomolecular 
domain is composed of the channel protein and the membrane and the 
solvent domain consists of extracellular (upper), channel pore (central), 
and intracellular (lower) regions. Particle species are indexed by 1, 2, and 3, 
for K+, CD, and H 2 O with radii ai = aK+, 0,2 = aci-; 03 = ohjO given in 
Table I. 

The SMIB method is an advanced method to treat singularities of protein 
charges and molecular surfaces [2^441145] . In SMIB, the electric potential gen¬ 
erated by the protein charges {qjS{r — rj) in ( 8 )) is modeled as a sum of an 
analytical Green function 0* in inhnite space and the Laplace potential 0° in 
the biomolecular domain with boundary values of 0* on dflm- The com¬ 
bined potential then dehnes an electric held — V((/>* + (fP) that acts on ions 
and water in the solvent domain from the molecular surface dVLm- The 
total potential 0 of all charged objects (ions, atomic charges, and polarized 
water) is then calculated by solving Eqs. ( 6 ) and (7) with the SMIB method 
for Eq. (7) across the interface dVLm of dielectric solvent and molecular 
domains. 
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TABLE I. Notations and Physical Constants 


Symbol 

Meaning 

Value 

Unit 

ks 

Boltzmann constant 

1.38 X 10-23 

J/K 

T 

temperature 

298.15 

K 

e 

proton charge 

1.602 X 10-^9 

G 

eo 

permittivity of vacuum 

8.85 X 10-1^ 

F/cm 

Cw 

water dielectric constant 

80 or 78.5 


Cm 

protein dielectric constant 

2 


e = e,(l-/2v2) 

dielectric operator, eg = e^eo 

in Eq. (6) 

F/cm 

e(r) ^ e 

dielectric function 

in Eqs. (6), (31) 


®Na+) ®Ca2+ 

particle radii 

0.95, 0.99 

A 

®K+) ®C1-) ®H20 

particle radii 

1.33, 1.81, 1.4 

A 

Ic 

correlation length 

1.2aK+ or 2 aQ^ 2 + 

A 

nB 

^K+ 

K+ diffusion coefficient 

1.96 X 10-5 

cm^/s 

nB 

^Na+ 

Na+ diffusion coefficient 

1.334 X 10-5 

cm^/s 

nB 

-^Ca2+ 

Ga^’*' diffusion coefficient 

0.792 X 10-5 

cm^/s 

^Cl- 

Gl~ diffusion coefficient 

2.032 X 10-5 

cm^/s 

nB 

H 2 O diffusion coefficient 

2.3 X 10-5 

cm^/s 


inside (outside) voltage 


V 


The molecular surface dVtrn as depicted in Fig. 1(a) is generated by rolling a 
probe ball (water molecule) with radius 1.4 A over a total of 554 spherical 
atoms in the GA protein [50]. In SMIB, the molecular surface is not hxed and 
is adaptively determined by the grid size so that the interface point is always 
in the middle of neighboring grid points. The resulting surface is thus free of 
geometric singularities. We refer to |15] for more details about the SMIB and 
surface generation methods. 

The NP equation (9) is then solved by the SG method for each particle species 
i once </> is known. An iterative process of solving PFl (6), PF2 (7), and 
NP equations is repeated again until convergent approximations of 0(r) and 
C'j(r) are found at all grid points. As noted above, convergence of this kind of 
iterative process is in general not guaranteed and must be checked at all grid 
points. We propose the following nonlinear iteration algorithm for the PNPF 
system (6), (7), and (9) using SMIB and SG methods: 
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Nonlinear Iteration Algorithm 1: 


1 . 


2 . 


Solve the Laplace equation —V^0‘^(r) = 0 in for the potential 0*’(r) 
once for all with 0 °(r) = 0 *(r) = |r — rjl) on dVLrn- 

Solve the Poisson equation —V • (eV0(r)) = 0 in ^2 at equilibrium for 
the initial potential with = 0 on dVL and the jump condition 


eV0 


Old 


■ n 


= + 0 °) ■ II 011 dVLrni where [u] denotes the jump 


3. 


unction across dll„^ ^^ 

Solve the PFl e, (/^V^ - 1 ) ^(r) = E,=i g*C'°'^(r) in 12, for 
y^New . n = 0 on dUm, = 0 on 512, and the Fermi distribution 

= In^^, F(r) = 1- 


pB 


4. 


-ed' 


New 


+ 


5. 


6 . 


C'Oi'i(r) = Cf exp (-A 0 Oi^(r) + ^*’'"(r)), 

Ef=i VjC^^'^ir). 

Solve the linearized PF2 —V ■ (eV0(r)) + p'(0°^'^)0(r) = 
p/( 0 Oid) 0 Oid equilibrium for the next potential 0 ^®’"(r) with the same 
jump and boundary conditions in Step 2 . Here p'(0) denotes the deriva¬ 
tive of the charge density functional p( 0 ) = exp (—/9i0 -1- 5*’’^) 

in 12 , with respect to 0 . 

Assign = a;pF0*^*‘^ + (1 ~ with a suitable relaxation param¬ 

eter cdpF and go to Step 3 if the error 0 ^®" — 0 ®^*^ in the inhnity norm 

CO 

is larger than a preset tolerance, else go to Step 6 . 

Solve the steady state NP equation — V-Jj(r) = 0 in 12, at nonequilibrium 
for Cj^®’"(r) and alH = 1, • • • , iP -|- 1 with Ji(r) = —Z2j [VC'i(r) -(- 0jC'j(r) 


V0Oid(r) _ C'^(r)V^*’'®fr) 


7. 

8 . 

9. 


^trc(j,) ^ In , C'j^®’"(r) = 0 on 512, and 

Jj(r) ■ n = 0 on 512^- 

Solve the PFl for v]>Ne’" as in Step 3 with in place of 

Solve the PF2 — V-(eV0(r)) = —(without linearization) at nonequi- 

librium for 0 ^®". 

Assign 0®^®^ = CdpNPF0^''^ A (1 — ^pnpf, 
parameter wpnpf and go to Step 6 if 
preset error tolerance, else stop. 


0 New suitable relaxation 

is larger than a 


/.New 


/.Old 


This is an SOR-like iteration algorithm modihed from that in [22]. The modi¬ 
fications include the additional solution processes at equilibrium in Steps 2 , 3, 
and 4, the extra PFl in Steps 3 and 7, Newton’s linearization for PF2 in Step 4, 
two relaxation parameters cjpf and (Upnpf in Step 5 and 9 with 0 < oipF, nipNPF 
< 1 (under relaxation), and the extra water NP equation —V ■ Jx+i(i’) = 0 in 
Step 6 . The stability and convergence rate are controlled by these two param¬ 
eters. If the parameter is close to zero, we will have more stable iteration but 
slower convergence. The correlation length 1^ and the Fermi distribution (or 
the steric potential S'*’’®) in Step 3 signify the difference between the classical 
PNP and advanced PNPF models. The stability and convergence are further 
complicated by these two physical properties for which a continuation method 
may be needed by introducing two stepping parameters Ac and Xs such that 
XJc and AgS'*’’® are gradually increased from Ac = A 5 = 0 toAc = A 5 = l [43] . 
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The water NP equation is not considered in previous PNP models and plays 
an essential role not only for numerical stability but also for physical validity 
because the void volume fraction P(r) = 1 — in Step 3 at any 

location r in the solvent domain hlg needs to be carefully checked in each itera¬ 
tion. Numerical errors in approximating the concentration functions Cj{r) for 
any particle species j = 1, ■■■, K + 1 could easily lead to an unphysical void 
fraction P(r) <0 at some r. Water molecules automatically adjust themselves 
in the PNPF model and move together with all ions in the system as the above 
iteration process converges to a stable and correct state, although the net water 
flow through the channel may be zero. Moreover, the water NP equation also 
dynamically determines the variable permittivity e(r)eo = es/(l + h(j')/p(j')) 
from the bath to the pore and thus automatically adjusts dielectric forces 
on ions along the channel pathway. These dielectric forces can have a decisive 
effect on biologically important conductance [HI] and on selectivity. For exam¬ 
ple, Na"*" vs. K"*" selectivity in Na"*" channels is only found when the dielectric 
function is handled in more detail |62ll63] . 

However, this SMIB-SG method and previous PNP methods suffer from a 
major difficulty in ion channel simulations. Those methods have difficulty in 
dealing with the essential property of selectivity, which of course is different in 
different types of channels with different structures. The L-type calcium chan¬ 
nel selects Ca^’*' over Na"*" of similar size and a potassium channel selects K"*" 
over Na"*" of the same charge. The following method is proposed to overcome 
this difficulty. 


B. A Multiscale Method for Calcium Channel 


Calcium channels have not yet been crystallized and so we use the Lipkind- 
Fozzard molecular model [ 51 ] of L-type calcium channels in which the EEEE 
locus (four glutamate side chains modeled by 8 0^/^“ ions) forms a high- 
affinity Ca^^ binding site that is essential to Ca^"*" selectivity, blockage, and 
permeation. Fig. 2(a) illustrates the binding site and the EEEE locus, where 
3 Ca^"*" are shown in violet, 8 in red, 2 H 2 O in white and red. Fig. 

2(b) is a cross section of a simplified 3D channel geometry for the present 
work, where the central circle denotes the binding site, the other four circles 
denote the side view of 8 0^^^“ ions, is the solvent domain consisting of 
two baths and the channel pore including the binding domain Deind, is 
the biomolecular domain with the boundary dflm, and dQ is the outside and 
inside bath boundary. Fig. 3 is a sketch of the binding site and 0^/^“ ions, 
where is the distance between the center of a binding Ca^’*' ion and the 
center Cj of any 0^/^“, and A is any point on the surface of the site. In our 
model, the 8 0^/^“ ions are not contained in the solvent domain Particle 
species are indexed by 1, 2, 3, and 4 for Na’*', Ca^+, CD, and H 2 O, respectively. 
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(a) 



(b) 


Fig. 2. (Color online.) (a) The Lipkind-Fozzard pore model, where 3 Ca^'*' are shown 
in violet, 8 0 ^/^“ in red, 2 H 2 O in white and red. Reprinted with permission from 
(G. M. Lipkind and H. A. Fozzard, Biochem. 40, 6786 (2001)). Copyright (2001) 
American Chemical Society, (b) A simplified Ca channel geometry in a cubic box 
with baths, pore, and binding site. The solvent region consists of two baths and 
the channel pore. The binding site fleind is contained in bnt the 0^/^“ ions are 
not in Qg. The outside and inside bath boundary is denoted by d^l. 


y 



Fig. 3. The binding distance between the center of the binding Ca^"*" ion and the 
center cj of the 0 ^/^“ ion is denoted by d.Q^' for j = 1 , • • • , 8 . A is any point on 
the surface of the binding ion. 


In [65] , we proposed an algebraic model for calculating the electrical potential 
(pfj and the steric potential in fiBind by using Coulomb’s law with the 
atomic structure of binding ion and atoms in a channel protein as shown in 
Fig. 3, without solving the Poisson-Fermi equation (8) in flBind- The volume 
of f^Bind is ^u unkuown variable Vf, that changes with different charges in the 
binding site. The algebraic model [65] dehned in f^Bind consists of the following 
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equations 


0\ = VbCf exp(-/3i0b + 5'^'') 
< = VbCf exp(-/324 + 

Ol = UftCf exp(S'!;''‘'), 

I. 



Vb - viOl - V2OI - U4O4 

- 


e 

ineo 



Zq1/2- 

Cj - ^1 


+ 


^l^Na+ 

®Na+ 


®Ca^+ j 


( 22 ) 


(23) 

(24) 


where Oj, O 2 , and O 4 denote the occupancy numbers of Na"*", Ca^"*", and H 2 O 
in Vb, respectively, 0 ^ and are average electrical and steric potentials, and 
\cj — A\ is the distance between A and Cj in Fig. 3. 


In this mean held, we allow 0\ and O 2 (and hence the total charge 0\ez-!^^+ + 
02 ezQg 2 +) to vary continuously subject to the condition on their sum O 4 +O 2 = 
1 in the binding volume Vb- Eqs. (22) and (23) uniquely determine the four 
unknowns Vb, O 4 , 0^ and with 0\ and O 2 being given. Eq. (24) uniquely 
determines the locations (cj) of 8 0 ^/^“ ions (and thus the binding distance 

or in Fig. 3) once 0^ is obtained. Note that the binding distance 
(or Cj) changes continuously with varying O 4 and O 2 but 0 ^ remains hxed, 
where the binding ion O^Na + O^Ca is a linear combination of Na"*" and Ca^’*'. 
Therefore, 0^^^“ ions are movable — the protein is hexible in our model — 
as their locations cj changes with varying O 4 and O^- 


For the half-blockage experimental condition [ 66 ] 

= Cf = 32 mM, ^^^ 2 + = Cf = 0.9 pM, (25) 

^ V 

Experimental Data 


we follow convention and assume relative occupancies of a hlled channel, 0\ = 
0.5 and O 2 = 0.5, and thereby obtain 0^ = —10.48 ksT/e, = —1.83, and 
Vb = 4.56 [65]. The binding experiments [ 66 ] used a hxed = Cf = 32 

mM and various Ca^^ bath concentrations <^^ 3 ^ 2 + = Cf that imply different 0\ 
and O 2 of Na"*" and Ca^’*' occupying the binding site. The occupancy numbers 
0\ and O 2 are determined by 

nb ^ _ nb _ nB 

^ - (92)0,)^, (26) 


where 0^ was just obtained from the case of equal occupancy. The occupancy 
ratio in (26) thus deviates from unity as Cf is varied along the horizontal axis 
of the binding curve from its midpoint value Cf = 0.9 /iM as shown in Fig. 5 
in [65] . 
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For nonequilibrium cases, the binding steric potential Sf, is assigned its equi¬ 
librium value in subsequent PNPF calculations, i.e., the void fraction F(r) in 
^Bind is assumed to remain unchanged from equilibrium to nonequilibrium. 
The electrical potential 0^ will be modihed by the membrane potential Id —14 
[13] and then used as a Dirichlet type condition for the potential function 0(r) 
in ffBind- For this multiscale method, the boundary conditions for the PF (8) 
and NP (9) equations are 

/ 

0 (r) = 0b(r) in fieind, d>(r) = 14,i on 
< C'i(r) = Cf on dQ, i = 1,2, 3,4, (27) 

Ji(r) ■ n = 0 on dflrn- 

.. 

Note that the electrostatic potential 0(r) is prescribed as a Dirichlet function 
0fe(r) whose spatial average in Dsind is the constant 0^. However, the binding 
domain Deind is treated as an interior domain instead of boundary domain for 
the NP equation (9). 

If a condition on the boundary is used to solve the Poisson (or PF) equation 
as in Algorithm 1, the resulting steric potential (as an output of 0(r) by 
(5)) may be incorrect in flBind because the atomic equations (23) and (24) 
are not used. We do not have any differential equation for the steric function 
^trc(r) fQj. which appropriate boundary conditions near Dsind can be imposed 
if a conventional method is used. The methods proposed in this paper are 
still coarse approximations to ion transport as the PNPF theory is in its early 
development. Nevertheless, the theory provides many atomic properties such 
as (23) and (24) that have been shown to be important for studying the binding 
mechanism in Cay channels [65] and are also important for the transport 
mechanism as shown in the next section. Incorporating atomic properties into 
continuum models is a step forward to improve and rehne the continuum 
theory. We refer to [451IB5] for more details of the algebraic model and its 
extension to PNPF. 

We summarize the PNPF solution process using the multiscale method as 
follows. 

Nonlinear Iteration Algorithm 2: 

1. Solve (23) and (24) for 0^ and in the binding site Dsind with the 
experimental data (25). 

2 . Choose any linear interpolation 0°*^ (an initial guess potential prohle) 
that links the binding potential 0 a to the zero potential at each bath 
boundary for the potential function 0(r). 

3. Solve the PB equation —V ■ (eV0(r)) = p(0°^'^) = at equilib- 
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4. 

5. 

6 . 

7. 


rium for with the Boltzmann distribution = Cf 


exp 




Compute the initial concentrations = C® exp . 

Solve the PFl e, - 1 ) ^(r) = in for 

the same conditions as in Algorithm 1 . 

Solve the linearized PF2 —V ■ (eV0(r)) + p'(0°^'^)0(r) = —+ 
in fls at nonequilibrium for 0'^®" with the conditions in (27). 
Solve the NP equation — V-Ji(r) = 0 in at nonequilibrium for C'f^®"(r) 


and alH = 1, • • • ,77 + 1 wit 
Go to Step 4 if 0^®" — 
error tolerance, else stop. 


1 the same conditions as in Algorithm 1. 


or 


^New _ (^Old 


is larger than a preset 


We do not need to solve the Poisson equation in the biomolecular domain 
VLrn that contains the singular charges of 8 0 ^/^“, since the effect of these 
charges on potentials has been included in the integral constraint we apply to 
the binding potential 0^ in (24). Consequently, we do not have to deal with 
the delta function in ( 8 ) and the potential jump conditions on dVLrn as used 
in Algorithm 1. The absence of jump conditions makes the approximation of 
PFl and PF2 more accurate since numerical methods for handling the jump 
conditions across molecular surfaces with singular cusps are subtle, complex, 
and thus prone to error [441145] . Moreover, the SOR-like scheme is not needed 
for this iteration. Application of the multiscale method to the NCX structure 
IST] is briefly discussed in [65] . It will be interesting to apply the method to the 
celebrated KcsA potassium channel [6H] and to recent structures of TRPVl 
[69] and CayAb channels P!- 

V. NUMERICAL RESULTS 


The main purpose of this work is to present numerical methods that are suit¬ 
able for continuum simulations of ion transport in different types of biological 
ion channels with particular interests in treating the excluded volume effect 
of all particles and the dynamical effect of water molecules. Numerical meth¬ 
ods are validated for accuracy with exact solutions of the PNP model for the 
GA channel. Numerical results of the PNPF model for both GA and calcium 
channels are all verihed with experimental data. 

A. Gramicidin A Channel 

The Scharfetter-Gummel method ( 21 ) is hrst validated with the following 
exact solutions for the PNP model [22] 


i COSTCOSUCOS^, r = {x,y, z) E VLrni 

K J m, ^ 28 ) 

cos X cos y cos z, r G 
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0 in Vtm, 


Gi(r) = 

(29) 

1 0.2cosxcos?/cosz -|- 0.3 in Vtg, 


0 in Dm, 


C^2(r) = 

(30) 

1 0.1 cosXcos 1 /cos -|- 0.3 in Vtg. 



Note that the right hand side of the Poisson equation in Algorithm 1 is not 
zero as the exact solution (28) has been imposed and the Green function 
0*(r) = (dTrCm |r — G'l) used in the jump condition on the 

molecular surface where the coordinates of the atoms in the GA 

channel protein are provided in the Protein Data Bank [52], the protein charge 
Qj and the radius of each atom j are obtained by the PDB2PQR software [7T] . 
and the total number of atoms is Q = 554. The optimal convergence (second) 
order, i.e. 0{h?), of the SMIB method for the nonlinear Poisson-Boltzmann 
equation has been conhrmed in |15]. The need for such validation has been 
pointed out before ua. It is easy to mistake convergence for accuracy in 
systems of PNP like equations [13] . 

For a full nonlinear PNP system (without steric, correlation, and water NP 
effects) using Algorithm 1, Table II shows that the optimal convergence order 
has been achieved for all PNP equations as well by using the SMIB method for 
the Poisson equation and the primitive FD method (12) for the NP equations 
in the nonlinear iteration process. In the table, errors are measured in the 
Loo norm. For example, 0.0927 = maxjjfc \(j){xi, yj, Zk) — where (pijk is the 
FD approximation of the Poisson equation and ^(xj, yj, Zk) is the exact value 
evaluated by (28) at the grid point (xj, yj, Zk) with the mesh size h = 1 A. The 
error tolerance for both linear solver and nonlinear iteration was set to 10 “®. 
All errors and orders (Ord) of convergence in Table II are similar to those in 
1221 , showing that the SMIB method in [15] is comparable to the original MIB 
1221 . 


When the primitive FD method is replaced by the SG method (21) for NP 
equations, it is surprising that the preset error tolerance 10 “® was satished by 
all SG approximations (pijk, and at all grid points for all mesh sizes 
as shown in Table III. Errors in Table III are much more smaller than those 
in Table II. This demonstrates that the SG is an optimal (exponential htting) 
method to discretize the NP equation as implied by the exact analysis of the 
ODE (20), since all solution functions in (28)-(30) are very smooth so that the 
assumptions made in (20) are valid. It only took 2 nonlinear iterations and 
about 1 hour and 8 minutes on a laptop computer with 2.6 GHz Intel GPU to 
reach the error tolerance for the case oih = 0.25 A. The corresponding matrix 
size is about 4.2 millions. The maximum potential difference A0j between 
any two adjacent grid points for the most coarse case (h = 1 A) is -1.045 
(not shown), which satishes the SG condition (18). This illustrates why the 
convergence has been achieved by the primitive FD without SG as shown in 
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Table II. 


TABLE II. Errors in 

norm 

by FD 


TABLE III. Errors by SG 

h{A) 

P 

Ord 

NPl 

Ord 

NP2 

Ord 

P 

NPl 

NP2 

1.00 

0.0927 


0.0505 


0.0211 


10 -® 

10 -® 

10-6 

0.50 

0.0245 

1.91 

0.0076 

2.73 

0.0042 

2.33 

10 -® 

10-6 

10-6 

0.25 

0.0060 

2.03 

0.0019 

2.00 

0.0010 

2.07 

10 -® 

10-6 

10-6 


We now study full PNPF (with steric, correlation, and water NP effects) sim¬ 
ulations of the GA channel using the SMIB and SG methods. Fig. 4 is a 
comparison of the I-V curves obtained by PNPF (lines) and the experimental 
data (symbols) from Gole et ah [TS] with bath K+ and Gl“ concentrations 
= 0.1, 0.2, 0.5, 1, and 2 M and membrane potentials AV = G — K = 0, 
50, 100, 150, and 200 mV. The PNPF currents in pico ampere (pA) were ob¬ 
tained with only one adjustable parameter, namely, the reduction parameter 6 
in the pore diffusion coefficients 6Df for all particle species, while all physical 
parameters in Table I were kept fixed throughout simulations. This kind of 
reduction parameter has been used in all previous PNP papers and is neces¬ 
sary in continuum simulations when compared with MD, BD, or experimental 
data because there is abundant qualitative evidence that the diffusion coef¬ 
ficient in channels is much smaller than in bulk, but quantitative estimates 
are not available, as well described by Gillespie in [SB] including Appendix 
and supporting material. In principle, all experimental data can be fitted by 
adjusting this parameter. For the PNPF currents at all G® and AV in Fig. 
4, we chose 6 = 1/4.7 which agrees with the range 1/3 to 1/10 obtained by 
many MD simulations of various channel models [2411741175] . indicating that 
the steric, correlation, and water NP properties have made PNPF simulations 
more closer (realistic) to MD simulations than previous PNP simulations for 
which the parameter 6 differs from MD values by an order to several orders 
of magnitude [24] . 

Furthermore, PNPF can also provide more physical properties that have not 
been observed by previous PNP models such as the variation of electric per¬ 
mittivity (dielectric function e(r) in Fig. 5) and water density (GH 2 o(r) in Fig. 
6 ) from bath to channel pore. Together with the electric (</>(r) in Fig. 7) and 
steric (S'*'''^(r) in Fig. 8) potentials, K’*' ions (in Fig. 9) are subject not only 
to the electric field — V0(r) but also to the steric (entropic) field VS'*'’‘^(r) as 
described in Eq. (10). These fields change with the variations of water density, 
other ion concentrations, voids F(r), and dielectric function e(r) at any loca¬ 
tion r in the solvent domain For example, the magnitude of electric fields 
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Fig. 4. A comparison of PNPF (lines) and experimental [73] (symbols) I-V results 
with bath K"*" and Cl“ concentrations = 0.1, 0.2, 0.5, 1, 2 M and membrane 
potentials AF = 0, 50, 100, 150, 200 mV. 

modified by the dielectric function e(r) can be as large as (80 —50)/80 = 37.5% 
of that by the constant permittivity 80eo for the bath condition = 2 M with 
the membrane potential 200 mV as shown in Fig. 5. The dielectric function 
in e(r)eo was calculated by 

e(r) =em + C'H 2 o(r)(ew - era)/C\Q (31) 

using the water density function CH 2 o(r) as proposed in [76]. The protein is 
most negatively charged around z = 13 A, where the pore is very narrow 
(about 1.6 A in radius) so that it is most crowded (most negative S'*’’^(r) = 
In in Fig. 8) there. The size effect of all particle species is clearly manifested 
by the steric function S'*'’‘^(r) in PNPF. These results provide one of the most 
comprehensive simulations on ion transport in real proteins using continuum 
models that we know of. 

The incompressibility of water and the mass conservation are important phys¬ 
ical properties that can be used to further verify continuum results. MD sim¬ 
ulations have shown that the GA channel can be occupied by two ions at 
moderately high concentration [77H75] . In Fig. 10(a), we observe that a total 
of 8 particles (water molecules plus ions) in the channel pore is conserved 
by PNPF but not by PNP as [KCl] increases from 0 to 2 M. The pore volume 
is determined by a length of 22 A (from 11 to 33 A in the channel axis in 
Fig. 9) and radii varying from 1.466 to 2.343 A along the axis (not shown). 
The PNPF water density profiles in Fig. 6 show that water molecules adjust 
self-consistently their configurations to accommodate ions (Fig. 9) in the 
two binding sites near the mouths of the channel as [KCl] increases. The com¬ 
plementary profiles of water and in Figs. 6 and 9 illustrate a continuum 
picture of six water molecules separating two ions in single file [75]. Figs. 
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Fig. 5. The averaged dielectric function ?(r) profiles at each cross section along the 
pore axis with C® = 0.1, 0.2, 0.5, 1, 2 M and AV = 200 mV. Figs. 5.3 - 5.6 are 
obtained with the same averaging method, C^, and AV. 



Fig. 6. The averaged water density CnjoCr) profiles. 

10 (a) and 10(b) also illustrate the saturation of ions and currents, respectively, 
as [KCl] increases. Note that PNP yields less ions (and hence currents) 
since the constant permittivity e = CwCo in Eq- (8 ) for PNP (with Ic = 0) is 
larger than the variable permittivity e(r)eo obtained by Eq. (31) for PNPF 
(with /(, 7 ^ 0) as shown in Fig. 5, i.e., larger e results in smaller charge density 
p (fewer ions) for the same 0. 

Therefore, the mass conservation and saturation results in Fig. 10 and the 
PNPF current results in Fig. 4 with the MD compatible parameter 9 = 1/4.7 
appear to justify the approximation formula (31) for calculating the variable 
e(r) that is an output for illustration. We emphasize that in our treatment, 
unlike most treatments of channels, dielectric and polarization effects are oper¬ 
ators that are outputs of the calculations. They are not assumed as constants. 


23 











Fig. 7. The averaged electric potential 4>{r) profiles. 



Fig. 8. The averaged steric potential profiles. 

The polarization effects of water are actually approximated by the dielectric 
operator e = 6^(1 — in Eq. (8) not by e(r). Obviously, the polarization 

effects (or equivalently the correlation effects of ions in PNPF) play a crucial 
role in very narrow channels that are more challenging to describe by classical 
PNP models or even by all-atom MD simulations as the current MD force 
helds do not include the electronic polarization effects [791180] . Of course, the 
single-hle picture by PNPF is still far from that by MD [78] due to inevitable 
averaging effects of numerous atoms in the system. Nevertheless, the electro¬ 
static, steric, and dielectric helds produced by PNPF may improve MD as well 
as continuum simulations in future studies. 


We make a hnal remark about the nonlinear iteration method for GA simu¬ 
lations. The two relaxation parameters of the SOR-like scheme in Algorithm 
1 were set to wpp = 0.3 and cjpnpf = 0.5 for all above results. The number 
of iterations for each PNPF I-V data point in Fig. 4 is given in Table IV. 
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Fig. 9. The averaged K"*" concentration Cj^+(r) profiles. 


(a) (b) 




Fig. 10. (a) Occupancy of H 2 O and K"*" in the GA channel pore by PNPF and PNP 
as [KCl] increases from 0 to 2 M. The total number of H 2 O and K'*' in the pore is 8 
m, which is conserved by PNPF but not by PNP (without steric and correlation 
effects), (b) Currents by PNPF and PNP at AV = 200 mV. PNP underestimates 
the currents. 


We do not need to solve NP equations when AV = 0. Iterations for solving 
PFl and PF2 in Steps 3 and 4 in Algorithm 1 increase with increasing bath 
concentrations as shown in the table. Iterations for solving K"*", Cl“, and H 2 O 
NP equations and then PFl and PF2 in Steps 6 , 7, 8 are all about 22 for 
G® = 0.1 to 1 M when AV 7 ^ 0. These numbers are more steady and less than 
those in [22] (see Table 7 in that paper). For = 2 M, iterations increase 
with increasing AV as those in [22]. Note that the relaxation parameter was 
set to different values for the first 3 steps in [22] whereas (Upp and (Upnpf were 
hxed throughout here. 
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TABLE IV. SOR Iterations 


C'^\AV 

OmV 

50 

100 

150 

200 

O.IM 

22 

22 

22 

22 

22 

0.2 

30 

22 

22 

22 

22 

0.5 

45 

21 

22 

22 

22 

1.0 

61 

21 

21 

21 

21 

2.0 

78 

34 

38 

43 

49 


B. Calcium Channel 


The calcium channel operates very delicately in physiological and experimental 
conditions as it shifts its exquisitely tuned conductance from Na’''-flow, to 
Na’'“-blockage, and to Ca^^-flow when bath Ca^’*' varies from trace to high 
concentrations. In |66], 19 extracellular solutions and 3 intracellular solutions 
were studied experimentally, in which the range of [Ca^’'']o is 10®-fold from 
to 10-2 M. 

PNPF results are in accord with the experimental data in [66] as shown in 
Fig. 11 under only the same salt conditions of NaCl and CaCl 2 in pure water, 
without considering other bulk salts in experimental solutions. With [Na’'']i = 
[Na’'']o = 32 mM and [Ca^+jj = 0, the membrane potential is fixed at —20 
mV (Vo = 0 and V = —20 mV) throughout, as assumed in Fig. 11 of [66] for 
all single-channel currents (in femto ampere fA) recorded in the experiment. 
The small circles in Fig. 11 denote the experimental currents from Fig. 11 of 
[66] and the plus signs denote the total currents calculated by PNPF, where 
the partial Ca^^ and Na"*" currents are denoted by the solid and dotted line, 
respectively. These two ionic currents show the anomalous fraction effect of 
the channel at nonequilibrium. The reduction parameter in 6Di was set to 
0 = 0.1 and all physical parameters in Table I were kept hxed throughout. 

Solution prohles of calcium channel are quite different from those of GA chan¬ 
nel as shown in Figs. 12 (dielectric function e(r)), 13 (water density GH 2 o(r)), 
14 (electric potential 0(r)), 15 (steric potential S'*'’'^(r)/cBT), and 16 (scaled 
flux densityl Jca^ 2 +(r)|). Fig. 12 displays the combined effects of correlation, 
polarization, and screening in this highly inhomogeneous electrolyte by means 
of the dielectric function e(r) in which water (Fig. 13) plays a more profound 
role than that of GA channel as the water density is dramatically reduced 
from 55.5 M in the bath to almost 0 in the binding site when [Ga^+jo = IQ-^ 
M. 
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Fig. 11. Single channel inward current in femto ampere (fA) plotted as a function 
of log^Q[Ca^'’']o. Experimental data of [BB] are marked by small circles and PNPF 
data are denoted by the plus sign and lines. 

Water is not allowed to occupy the binding site because Ca^^ occupies it in 
this bath condition and the 8 ions in the EEEE locus are electrically 

attracted toward the binding Ca^’*' as illustrated by Fig. 2(a). The EEEE locus 
is very hydrophobic in this condition. Without using the atomic properties of 
water and ions inside the solvent domain like those in (22)-(24), continuum 
models are not likely to produce results like Figs. 12 and 13. Mathematically, 
the Dirichlet condition 0(r) = 0;,(r) in the interior of his, he. fleind C in 
(27), is crucially important to connect the continuum Poisson-Fermi model (4) 
in r2s\flBind to the molecular (Coulomb) model (22)-(24) in fleind- From the 
binding formula (24), the pore radius is enlarged by the binding Na"*" when 
[Ca^’'']o decreases from 10“^ to 10“^'^ M, i.e., Na^ occupancy (Oi in (26)) 
increases. The enlarged radius allows more space for water molecules in the 
channel pore as shown in Fig. 13. 

Novel Features. The steric potential profiles shown in Fig. 15 represent 
the novelty of the PNPF theory. All effects of volume exclusion, interstitial 
void, configuration entropy, short range interactions, correlation, polarization, 
screening, and dielectric response of this nonideal system are described by the 
steric functional S'*'’‘^(r). The steric potential in the binding region decreases 
drastically from —1.30 to —10.34 ksT as [Ca^''']o increases from 10“^'^ to 10“^ 
M. However, the electric potential remains almost unchanged as shown in Fig. 
14 following the linear occupancy model (26). In physiological bath conditions 
[Ca^’'']o = 10“^ ~ 10“^ M, Fig. 13 shows that the region containing the bind¬ 
ing site with the length about 10 A is very dry (hydrophobic), which agrees 
with the recent crystallographic analysis |S7| of the Ca^’*' binding site of the 
related protein NCX with the EETT locus showing a hydrophobic patch (also 
about 10 A in length) formed by the conserved Pro residues. The hydropho- 
bicity near the binding site in our model is described by the continuous water 
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Fig. 12. The averaged dielectric function ?(r) profiles at each cross section along the 
pore axis for various [Ca^'’']o ranging from 10“'^'^ M to 10“^ M. All the following 
figures are obtained with the same averaging method and the same range of [Ca^'’“]o. 



Fig. 13. The averaged water density C'h 2 o(i') profiles. 

density function C'H 2 o(r) via the continuous steric function namely, 

the Fermi distribution Ch 2 o(i') = C'hjO in (5). At [Ca^’'']o = 10“^ 

M, the magnitude of the steric energy = —10.34 is comparable to 
that of the electrostatic energy 0 = —10.48 ksT/e. This surprisingly large 
energy due only to the steric effect has not been quantihed and observed by 
other continuum methods in Cay channel modeling, as far as we know. 

The variable steric potential with respect to bath concentrations as shown in 
Fig. 15 plays a similar role as the conhnement potential in MC simulations 
for constraining the 8 mobile 0^/^“ ions of protein glutamates in a hlter |59j . 
These are two different approaches to modeling the flexible glutamates and the 
steric effect of ions. The excluded volumes of electrolyte and glutamate ions 
are explicitly given as an input in MC simulations by using the conhnement 
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Fig. 14. The averaged electric potential 4>{r) profiles. 



Fig. 15. The averaged steric potential profiles. 

potential in a fixed filter whereas the volumes are implicitly calculated in 
PNPF simulations and are outputs that describe the steric potential in a 
variable binding site. We had difficulties obtaining nonequilibrium results as 
those in Fig. 11 in our early attempt to use a hxed conhnement potential in 
a hxed hlter partly because the conhnement potential would generate large 
artihcial electric helds near the boundary of the hlter in continuum setting. It 
is also difficult to incorporate the conhnement potential into the hux density in 
Eq. (10) since the conhnement potential is hxed and cannot not be explicitly 
decomposed to the electrostatic and non-electrostatic parts the way and 
do in (10). These difficulties are typical of inconsistent calculations. Imposing 
potentials (whether in simulations or theories) requires injection of charge 
and mass that is not present in the real system. The injection occurs at a 
sensitive part of the system, the selectivity hlter. The approach here avoids 
these difficulties. 
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Fig. 16. The averaged flux density |J(j^ 2 +(r)| profiles. 


As observed from Fig. 11, ionic transport is blocked by the competition be¬ 
tween Na"*" and Ca^"*" ions in the range [Ca^’'']o = 10“^'^ ~ 10“"^'^ M. In this 
blocking range, the corresponding steric prohles in Fig. 15 are wider indicating 
that the water density or the void volume is more evenly distributed. From 
Fig. 11, we observe that Ca^’*' currents increase dramatically when [Ca^’'']o 
increases from 10“^'^ to 10“^ M in the physiological mM range of the channel, 
so does the corresponding flux density as shown in Fig. 16. 

Numerical Verification. All the above results were obtained by using the 
standard FD method (see 03 for instance) for the Poisson-Fermi equations 
(6)-(7) and the Scharfetter-Gummel method (21) for the flux equation (9). We 
now provide more numerical details for the extended SG stability condition 
(18) and explain why the primitive method (12) fails. We hrst verify the SG 
method at equilibrium (AG = G = id = 0) for which the PNPF solution 
should agree with the PF solution as shown in Table V, where the averaged 
Na’*' and Ga^+ concentrations in the Alter (binding site) are denoted by Gf 
and Gj, respectively. Here, the approximate PF solution of Gj(r) in (5) is 
treated as an exact solution to justify the approximate PNPF solution. The 
PNPF concentrations agree quite well with the PF concentrations indicating 
that the SG method (21) works well. Note that Gf and Gf are all bounded by 
their respect maximum values 462.39 and 408.57 M at very large electrostatic 
potential 0 = —10.48 in the Alter, as guaranteed by the Fermi distribution (5). 
We use the SOR method [8T] for solving all linear algebraic systems with the 
relaxation parameter taken to be 1.9. The error tolerance of the SOR linear 
solver is 10“® because the boundary bath condition [Ga^’'']o for (21) is in the 
10®-fold range. The error tolerance for solving each PDF in the PNPF model 
in Algorithm 2 is 10“^. 
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TABLE V. Verification of the SG method (21) 


[Ga2+] in M 

0.9- 10-® 

10“^ 

10-^ 

10-3 

10-2 

PF Cf/C^ 

PNPF 

61.7/63.7 

61.7/63.4 

10.1/116.9 

10.1/116.1 

1.1/126.2 

1.1/126.2 

0.11/127.3 

0.11/127.3 

0.01/127.4 

0.01/127.4 


Primitive FD method fails. We next look more closely into the numerics 
of the SG discretization concerning the SG condition (18) at nonequilibrium 
(AV = —20 mV) under the conditions as those in Fig. 11. In Table VI, —/SjA^ 
denotes the maximum difference of —/9j0(r) between all adjacent pairs of grid 
points in 3D, where the subscript i denotes the ionic species not the grid node. 


Since 


Ga 


2 + 


varies in the 10®-fold range, the maximum difference —/SjA^ 


varies in a range for each ionic species i as shown in the table. Moreover, the 
adjacent pair of grid points at which the maximum difference is obtained may 
differ for different NP equations. The other two maximum differences AS'*'''^ 
and —/SjA^ + AS'*'''^ are similar dehned. Note also that the three maximum 
differences may occur at different pairs of adjacent grid points with the mesh 
size h = 1 A even for the same NP equation. 


TABLE VI. Ghecking the SG condition (18) 



-l3iA(j) 

A^trc 

-AA0 + 

NPl (Na+) 

2.43 ~ 2.79 

0.34 ~ 1.75 

1.16 ~ 2.51 

NP2 (Ga2+) 

4.86 ~ 5.59 

0.34 ~ 1.75 

4.03 ~ 5.30 

NP3 (G1-) 

2.47 ~ 2.82 

0.34 ~ 1.75 

0.34 ~ 1.75 


From Table VI, we observe that the primitive FD method (12) violates the 
SG condition (18) for all NPl, NP2, and NP3 (without and for NPl and 
NP2 (with S'*’’^). The worst case of the violation occurs in the Ga^"*" flux equa¬ 
tion NP2 with or without as analyzed in (17). Obviously, the primitive 
FD is not suitable for PNPF simulations on Gay channels. The SG not only 
delivers stable results for all equations in the PNPF model at all experimental 
conditions but also converges very rapidly in the nonlinear iteration. 

VI. CONCLUSION 

The classical Scharfetter-Gummel (SG) method for semiconductor devices was 
extended to include the steric potential for biological ion channels in this pa¬ 
per. The steric potential — a key feature of the PNPF theory — represents 
a combination of various effects of volume exclusion, interstitial void, config¬ 
uration entropy, short range interactions, correlation, polarization, screening. 
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and dielectric response in a complex fluid system of ion channel. The simpli¬ 
fied matched interface and boundary (SMIB) method together with the SG 
method was shown to yield optimal convergence for a PNP model with exact 
solutions of the gramicidin A channel. The primitive finite difference method 
without SG was shown to lead to unphysical approximations for an L-type cal¬ 
cium channel due to the violation of the generalized SG condition presented 
here. Two algorithms based on the SMIB and multiscale methods have been 
presented for these two different types of channels depending on whether water 
is allowed to pass through the channel pore. Numerical results for both chan¬ 
nels are in accord with the respective experimental results. Gompared with 
previous PNP models, new physical details by PNPF such as water dynamics, 
dielectric function, voids, and steric energy in the system have been illustrated 
and discussed. The PNPF model differs from most channel models in several 
respects. It computes dielectric properties as an output that in fact vary with 
position and with experimental condition. It provides a fourth order partial 
differential equation to describe current flow, of the general Gahn-Hillard type, 
which has a richness of behavior beyond the usual second order PNP descrip¬ 
tion. Practically, it is important that PNPF is much easier to compute in three 
dimensions than other steric PNP models. 
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